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Abstract 



Using molecular dynamics computer simulations we investigate how the glass 
transition and the properties of the resulting glass depend on the cooling rate 
with which the sample has been quenched. This is done by studying a two 
component Lennard-Jones system which is coupled to a heat bath whose tem- 
perature is decreased from a high temperature, where the system is a liquid, 
to zero temperature, where the system is a glass. The temperature of the 
heat bath is decreased linearly in time, i.e. T;, = Tq — 7t, where 7 is the 
cooling rate. In accordance with simple theoretical arguments and with ex- 
perimental observations we find that the glass transition, as observed in the 
specific heat and the thermal expansion coefficient, becomes sharper when 7 
is decreased. A decrease of the cooling rate also leads to a decrease of the 
glass transition temperature Tg and we show that the dependence of Tg on 
7 can be rationalized by assuming that the temperature dependence of the 
relaxation times of the system is given by either a Vogel-Fulcher law or a 
power-law. By investigating the structural properties of the glass, such as 
the radial distribution functions, the coordination numbers and the angles 
between three neighbor-sharing particles, we show how the local order of the 
glass increases with decreasing cooling rate. The enthalpy and the density of 
the glass decrease and increase, respectively, with decreasing 7. By investi- 
gating the 7 dependence of clusters of nearest neighbors, we show how these 
observations can be understood from a microscopic point of view. We also 
show that the spectrum of the glass, as computed from the dynamical matrix, 
shows a shift towards higher frequencies when 7 is decreased. All these eff'ects 
show that there is a significant dependence of the properties of glasses on the 
cooling rate with which the glass is produced. 

PACS numbers: 61.43.Fs, 61.20.Ja, 02.70.Ns, 64.70.Pf 
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I. INTRODUCTION 



If a liquid is cooled rapidly enough, so that the crystallization at or slightly below the 
melting point is avoided, the final state of the material will be an amorphous solid, i.e. a 
glass. This so-called glass transition, i.e. the transition from a liquid state to an amorphous 
solid state, is essentially the falling out of equilibrium of the system, because the relaxation 
times of the system at the glass transition temperature exceed the time scale of the experi- 
ment. Therefore the resulting glass is not in thermal equilibrium. Thus it can be expected 
that the properties of the glass will in general depend on the production history of the glass, 
as, e.g., the cooling rate with which the sample was cooled or, if the glass transition is pres- 
sure driven, on the rate with which pressure was applied. That such dependencies indeed 
exist has been shown in various experiments in which it was investigated how the density 
of the glass, the glass transition temperature, the specific heat or the Mossbauer spectrum 
depend on the cooling rate |I|J^. 

The dependence of the properties of a structural glass on its production history have 
also been studied in computer simulations. Fox and Andersen have investigated for a one- 
component Lennard- Jones system the dependence of the density on the cooling rate and 
Baschnagel et al have studied on how various properties of a polymer glass depend on the 
cooling rate ||^. However, because of their relatively complex model, the authors of this last 
work were only able to vary the cooling rate by two decades. Lai and Lin investigated the 
dependence of the Wendt-Abraham parameter on the cooling rate and Miyagawa and 
Hiwatari gave evidence that for their binary soft sphere model the structural properties of 
the glass are independent of the cooling rate, that, however, dynamical quantities depend 
on the cooling rate 0. Also Speedy investigated this type of question and reports that the 
rate with which a system of hard spheres is compressed does not affect the density of the 
glass 0]. 

Thus although computer simulations are in principle a very useful tool in order to in- 
vestigate the production history dependence of the properties of glasses, since they allow 
to study a glass in its full microscopic details, relatively few simulations have been done to 
address this question in detail. One reason for this is that the effects to be expected are 
relatively small, usually on the order of a few percent when the cooling rate is varied by 
two or three decades []T],0], and that therefore the statistical accuracy of the data has to be 
rather high in order to detect these effects at all. However, as the cited papers have shown, 
this problem is not insurmountable and since the potential payoff of such types of computer 
simulations, namely to understand on a microscopic level how the production history of the 
glass affects its macroscopic properties, seems to justify the effort involved. 

A further reason for investigating the dependence of the properties of the glass on the 
cooling rate originates from those computer simulations in which a relatively realistic type of 
model is used in order to simulate real glass formers, such as Si02. In these simulations the 
quahty of the applied model is in most cases judged by its ability to reproduce some exper- 
imentally determined quantities, such as the temperature dependence of the density, radial 
distribution functions or bond angle distributions functions 0. Before such a judgement 
can be made, however, it is necessary to see how the properties of the glass, as determined 
from a computer simulation which involves cooling rates that can be more than ten orders 
of magnitude larger than the ones in real laboratory experiments, are affected by these high 
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cooling rates. It might well be, that a discrepancy between the prediction of a model and 
the experimental data is not due to the inadequacy of the model, but due to the too high 
cooling rate with which the glass on the computer was produced. Thus it is important to 
separate these two possible sources of the discrepancy, and investigating the cooling rate 
dependence of the properties of the glass is one possible way to do this. 

In the present paper we investigate the dependence of the glass transition phenomenon 
and the properties of the resulting glass on the cooling rate. The system we use is a binary 
mixture of Lennard- Jones particles, since this is a prototype of a simple glass former. The 
goal of this work is to show that for this system there are significant dependencies of the 
details of the glass transition and of the properties of the glass on the cooling rate. Fur- 
thermore we try to understand these dependencies from a microscopic point of view, thus 
making use of one of the advantages of computer simulations. 

The rest of the paper is organized as follows. In Sec. II we introduce our model and 
give some of the details of the simulation. In the first part of Sec. Ill we present our results 
regarding the cooling rate dependence of the glass transition and in the second part of that 
section we give the results of the coohng rate dependence of the resulting glass. In Sec. IV 
we then summarize and discuss these results. 



II. MODEL AND DETAILS OF THE SIMULATION 

The system we are investigating is a binary mixture of particles, subsequently called 
type A and type B, each of them having the same mass m. The interaction between two 
particles of type a and (3 {a, P E {A,B}) is given by a Lennard- Jones potential, i.e. V^/3(r) = 
4ea/3[(?"/fa/3)^^ — / cTa/sY]. The parameters ea/3 and aa/3 were chosen to have the following 
values: aAA = 1-0, caa = 1-0, <jab = 0.8, e^s = 1.5, asB = 0.88 and ess = 0.5. In order to 
decrease the computational burden this potential was truncated and shifted at r = 2.5aai3- 

The model presented has already been used to investigate the dynamics of strongly 
supercooled liquids ||^. From these previous investigations it is known that this system is a 
good glass former, in that it is not prone to crystallization, even if the system is cooled very 
slowly. Therefore this system is a good model for investigating the dependencies of various 
properties of the glass on the cooling rate. 

In the following all results are given in reduced units. The unit of length is aAA, the unit 
of energy is caa and the unit of time is {ma\j^/ ASeaaY^'^ ■ In order to allow the investigation 
of the dependence of the density of the glass on the cooling rate and the computation of the 
specific heat at constant pressure, the simulations were performed at constant pressure Pext- 
This was done by using the algorithm proposed by Andersen |T^. The external pressure 



Pext was chosen to be 1.0 and we used the value of 0.05 for M, the mass of the piston. 

In order to simulate the cooling process we proceeded as follows. The system was equi- 
librated in a iVpT-ensemble at a high temperature Tq which was for most cooling processes 
chosen to be 2.0 (see below for a modification of this value). At this temperature the relax- 
ation times are short and thus equilibration is fast. After this equilibration period we began 
with the cooling process. For this we integrated the equations of motion with the velocity 
form of the Verlet algorithm, using a time step of 0.02. Every 150 time steps the velocities 
of all the particles were replaced by new velocities which were drawn from a Boltzmann 
distribution corresponding to the temperature Tf, of the heat bath. The temperature of this 
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bath was decreased linearly in time, i.e. Tb{t) = Tq — 7^, where 7 is the cooling rate. This 
cooling process was continued until the temperature of the bath was zero, i.e. for a time 
T0/7, which means that for the smallest cooling rate (7 = 3.125 ■ 10~^) the length of the run 
was 8 ■ 10^ time steps. Using the multi conjugate gradient method the coordinates of 
the particles as well as the volume were subsequently allowed to relax to the nearest local 
minimum of the potential energy hypersurface. It should be mentioned that an equivalent 
way to relax the configurations would have been to continue the molecular dynamics simu- 
lation at Tfe = for a long time and that the only reason for using the conjugate gradient 
method was computational efficiency. 

The cooling rates we investigated were: 7 = 0.02, 0.01, 0.005, 0.0025, 0.001, 5.0 ■ 10-^ 
2.5 ■ 10-^ 1.0 ■ 10-^ 5.0 ■ 10~^ 2.5 • 10~^ 1.25 ■ 10-^ 6.25 ■ lO-^ and 3.125 ■ 10"^ The number 
of A particles was 800 and the one for the B particles was 200. In order to test for finite 
size effects we also did some runs for a system twice as large. In order to decrease the error 
bars of our results we averaged the results for each value of 7 over ten independent runs, 
each of which was equilibrated at T = 2.0 for 2000 time units. Since at the smallest cooling 
rate each of these runs took about 280 hours of cpu time on a IBM RS6000/375 it is today 
probably not possible to investigate a range of cooling rates that is significantly larger than 
the one investigated here. 

Finally we make a remark on the choice of the starting temperature Tq for the cooling 
process. It can be expected that for a given cooling rate 7 the temperature Tg{'~f) at which 
the system falls out of equilibrium, and hence also the properties of the resulting glass, is 
independent of the starting temperature Tq, as long as Tq is sufficiently above Tg, or, to 
put it differently, that the time constant of the cooling, [{—dT/dt)/To]~^ = I/7T0, is larger 
than the intrinsic relaxation time of the system at the temperature Tq. In order to save 
computer time we therefore chose for the small values of 7 a value for Tq which was lower 
than the corresponding value for large values of 7. In particular we chose Tq to be 2.0 for 
7 > 1.0 ■ 10-^ equal to 1.0 for 7 = 5.0 ■ 10-^ 2.5 ■ 10-^ 1.25 ■ 10"^ and equal to 0.5 for the 
two smallest values of 7. In all cases we made sure that the value of the glass transition 
temperature is indeed smaller than Tq, which was done by investigating the temperature 
dependence of the enthalpy and the density (see below for details). 



More information regarding the details of the simulation can be found in Ref. [12 



III. RESULTS 

In this section we will report the results of our simulation. In the first subsection we will 
discuss the coohng rate dependence of the system during the cooling process, i.e. look at 
the system at finite temperatures, and hence investigate the details of the glass transition. 
In the second subsection we will present the results on the cooling rate dependence of the 
so obtained glass at T = 0. 



A. Finite Temperatures 

In this subsection we will investigate the properties of the system during the cooling pro- 
cess and see how these properties depend on the cooling rate. One of the simplest quantities 
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one can investigate is the enthalpy H(Tb) = Ekm + Epot + PextV + MV'^/2, where V is the 
volume of the system. This quantity has been studied in many previous computer simula- 
tions and it was found that on lowering the temperature a noticeable bend in H{Th) occurs 
at the temperature at which the systems falls out of equilibrium. Therefore the temperature 
at which this bend occurs is usually identified with the glass transition temperature Tg. 

In Fig. we show the enthalpy if as a function of T^, the temperature of the bath, for 
all cooling rates investigated. We see (Fig. |^a) that at high temperatures the curves for the 
intermediate and small cooling rates (bottom curves) fall onto a master curve. This master 
curve is given by the equilibrium value of the enthalpy of the system in its (supercooled) 
liquid-like state, i.e., by the curve one would obtain for an infinitely small cooling rate. 
For a given cooling rate the enthalpy of the system will be given by this equilibrium curve 
provided that the temperature of the system is high enough to allow it to equilibrate before 
a significant change in temperature, due to the cooling, has happened. If this condition is 
violated the system will fall out of equilibrium, and start to show a solid-like behavior, which 
is reflected by the bending of the curve for the enthalpy. Thus the system has undergone 
a glass transition. We therefore expect that the temperature Tg at which this transition 
happens is smaller the lower the cooling rate is. That this is indeed the case is shown in 
Fig. |I]b in which we show the temperature region in which the glass transition occurs on an 
enlarged scale. We recognize that the smaller the cooling rate is (curves at the bottom) the 
later these curves split off from the liquidus curve. Later we will investigate this phenomenon 
more quantitatively. 

From Fig. |l]a we also recognize that for the fastest cooling rates (top curves) the system 
falls out of equilibrium already at the starting temperature Tq. These cooling rates are so 
large, that it is no longer appropriate to see the cooling as a process which ultimately forces 
the system to fall out of equilibrium because its relaxation times exceed the time scale of the 
simulation, but rather as similar to a steepest descent procedure of the system with respect 
to the enthalpy. Therefore we can expect that for such large cooling rates the temperature 
dependence of various quantities is qualitatively different from the one for smaller cooling 
rates and the fact that the curves for the enthalpy split off from the master curve already 
at Tq is an example for such a different type of behavior. 

The just presented behavior for the enthalpy is very similar to the one found for the 
density p. The temperature dependence of p (see Fig. |^) shows that p{Tb), similar to H(Th), 
has a bend in the vicinity of ^ 0.4 and that the curves corresponding to the various 
values of 7 split off from the liquidus curve at a temperature that is lower the smaller 
7 is. We also note that for this system the temperature dependence of the density is a 
monotonous function of temperature. This will allow us later to define in a simple way a 
glass transition temperature. It has to be remarked that this is not a universal property, 
since other systems, such as Si02, show an anomaly in the density and thus make the 



extraction of a glass transition temperature from the density much more difficult We 
also recognize from the inset in Fig. ^, that for this system the density varies by more than 
a factor of two when the temperature is changed from T = 2.0 to T = 0, which is also in 
strong contrast with the case of the network glass former Si02 ||T3| . 

By differentiating the enthalpy with respect to the temperature we obtain Cp, the specific 
heat at constant pressure. Since our data for H(Th) was a bit too noisy to compute its 
derivative directly from a difference quotient, we approximated H(Th) by a sphne under 
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tension [jl4[ and then determined Cp by taking the difference quotient of this sphne. In Fig. |^ 
we show the so obtained curves for Cp for all cooling rates investigated. We see (Fig. ^ja) 
that at high temperatures and very fast cooling rates Cp drops from a value around 5 to 
significantly smaller values when the temperature is increased. This decrease is again a 
signature that for these large cooling rates the system falls out of equilibrium immediately 
after the cooling process is started. Thus this decrease should not be seen as the generic 
behavior of Cp for intermediate and small cooling rates. Rather we see that for such values 
of 7, Cp is almost independent of Ti, for larger than 0.8. Only when the temperature 
is decreased to approximately 0.5 we see that the value of Cp starts to drop from a value 
around 5.0 to a value around 3.0, the value expected for a classical harmonic system. We 
also recognize that even at these low temperatures the values for Cp for the fast cooling rates 
are considerably larger than the harmonic value 3.0. This shows that for these fast quenches 
the system behaves very anharmonic even at very low temperatures. 

In order to see the dependence of Cp on the cooling rate clearer we show in Fig. |^ the 
temperature range in which Cp shows the strong decrease on an expanded scale. We recognize 
that the temperature range in which this drop occurs is identical to the temperature range 
in which the bend is observed in the enthalpy and the density (see Figs. |l| and ||) and thus 
is related to the glass transition. 

From Fig. |^ we recognize that the drop of the curves becomes sharper the smaller 
the cooling rate is. Thus we find that the glass transition becomes more pronounced with 
decreasing cooling rate, in accordance with the simple argument put forward by Angell et al 



T5| which says that, since the relaxation times are a strong function of temperature, e.g., of 



Arrhenius or Vogel-Fulcher type, the temperature difference T2 — Ti for which the relaxation 
times change by, say, a factor of 100, i.e. r(Ti)/r(T2)=100, decreases with decreasing Ti. 
Since we have seen that the glass transition temperature Tg decreases with decreasing 7, 
it thus follows that the transition happens in a temperature range Ti < Tg < T2 which 
becomes narrower the smaller the cooling rate is. 

A qualitatively similar behavior to the cooling rate dependence of Cp is found for the 
cooling rate dependence of the thermal expansion coefficient ap{Th). This quantity is given 
by ap{Tb) = dV(Tb)/dTh\^ and is shown in Fig. ^ for all cooling rates investigated. 
To compute the derivative we made again use of splines under tension. The drop in Op 
at high temperatures for large values of 7 should again be viewed as the atypical behavior 
of ctp for very large cooling rates. For intermediate and small cooling rates we see that 
Op shows an almost linear dependence on temperature until it shows a rather sudden drop 
when the glass transition temperature is reached. As in the case of the specific heat this 
drop becomes sharper with decreasing cooling rate. Thus we see also for this quantity that 
the glass transition becomes more pronounced the smaller the cooling rate is. However, it is 
clear that even for our smallest cooling rate the glass transition is still smeared out over a 
relatively large temperature range, unlike the real experiment at much smaller cooling rates. 

From Figs. ^ and ^ we see that the temperature at which the curve for a given value 
of the cooling rate splits off from the liquid-like branch decreases with decreasing 7. This 
means that the glass transition temperature Tg decreases with decreasing cooling rate. In 
order to investigate this effect in more quantitative terms we use the concept of the "fictive 



temperature" as introduced long ago by Tool and Eichlin [Id]. This concept is founded 



on two observations: The first is that at high temperatures the temperature dependence 
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of most bulk properties of the system is independent of the coohng rate, i.e. the curves 
for the different coohng rates fall onto a master curve, as we see it, e.g., in Figs. |l| and 
0. This master curve can thus be approximated by a low order polynomial and therefore 
extrapolated to low temperatures. The second observation is, that the low temperature 
behavior of these bulk properties, is, apart from an additive constant, independent of the 
cooling rate (see, e.g.. Figs. and |^). The physical reason for this is that the temperature 
dependence of such bulk quantities is not very sensitive on the details of the microscopic 
structure of the glass, which themselves are dependent on the cooling rate. Therefore also 
this low temperature branch can be approximated by a low order polynomial and thus 
extrapolated to higher temperatures. The point where the extrapolations of the curves at 
high and low temperatures intersect is the fictive temperature which in this work we will 
call the glass transition temperature Tg. Note that this construction is only possible if the 
curves actually fall onto a master curve at high temperatures. Since this is not the case for 
the fastest cooling rates investigated here, we determined Tg only for the intermediate and 
small cooling rates, i.e. 7 < 0.005. 

In order to determine the low temperature master curve we shifted the curves for the 
different values of 7 vertically by a 7 dependent amount until the curves lay on a master 
curve in the temperature range < < 0.2. The so obtained master curve was fitted in 
this temperature range with a quadratic polynomial. (We note that a linear function, as 
has been used in several other investigations, did not give a good fit to our low temperature 
data.) The master curve for the high temperature data in the range 0.6 < Tb < 0.8 was 
fitted well by a linear function in the case of the enthalpy, whereas a quadratic polynomial 
was needed in the case of the density. 

In Fig. ^ we present the glass transition temperatures, as determined by the procedure 
described above, as a function of the cooling rate. The error bars are estimated from the 
uncertainty in the determination of the intersection points of the extrapolated curves. We 
see that, as expected, Tg decreases with decreasing cooling rate. We also recognize that 
the values of Tg as determined from the enthalpy are not exactly the same as the ones 
as determined from the density, in that, e.g., for not too large cooling rate, the former 
are systematically lower than the latter. However, at the moment we are not able to say 
whether this behavior has some real underlying physical reason or whether it just refiects 
some systematic error in our determination of Tg. It is likely that this difference becomes 
smaller at very small cooling rates, where the glass transition is much sharper than in the 
present case. 

The dependence of Tg on the cooling rate has also been determined in experiments 0. 
It was found that the functional form 

is able to give a satisfactory fit to the data when 7 is varied over 3 decades. The functional 
form of Eq. (P can be rationalized by assuming that the relaxation time r(T) of the system 
shows a Vogel-Fulcher dependence on temperature, i.e. r(T) = Aexp{B /(T — Tq)) and 
by assuming that the system falls out of equilibrium at that temperature Tglj) where the 
relaxation time is of the order of the inverse of the cooling rate, i.e. r(Tg(7)) oc 7"^. 
In Ref. it was shown that, at constant volume, the temperature dependence of the 
relaxation times of the system is also fitted very well by a power-law, a functional form 
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proposed by the so-called mode-coupling theory |T^. If we use this functional form for r, 
i.e., t(T) = A/{T — TcY, we expect that Tg depends on 7 via 



T,(7) = T, + (^7)^/^ . (2) 

Also included in Fig. | are fits to the data with the functional forms given by Eq. (|I]) (solid 
curves) and Eq. (0) (dashed curves). We see that both fits are able to represent the data very 
well. The temperature Tq from Eq. (|I|), i.e., the glass transition temperature for an infinitely 
slow cooling process, are 0.334 and 0.348 for the enthalpy and the density, respectively and 
the corresponding values for T^. from Eq. (|^) are 0.378 and 0.386, respectively. Also in this 
case it is not clear whether the difference between the values of Tq (or Tc) for H and p 
has some underlying physical reason or whether it is just a systematic error in the way we 
determined Tg{j). 

We also note that the estimated value for is significantly larger than the one estimated 
for Tq. The reason for this is that the critical temperature of mode-coupling theory is 
expected to overestimate the glass transition temperature since the theory neglects certain 
types of relaxation processes which restore ergodicity. Thus in experiments it is found that 
at the temperature T^. the system shows only a smeared out singularity and that Tc> Tq. 



B. Zero Temperatures 

After having discussed in the preceding subsection the cooling rate dependence of various 
quantities at finite temperatures we now investigate how the properties of the glass at zero 
temperature depend on the cooling rate. 

We have seen in Fig. |l] that the smaller the cooling rate is the lower is the temperature at 
which the curve for the corresponding 7 splits off from the liquidus curve. Thus we expect 
that the value of the enthalpy at zero temperature will depend on 7. That this is indeed 
the case is shown in Fig. where we plot if/, the final enthalpy after the quench (and the 
subsequent relaxation of the system, see Sec. |1|) as a function of 7. As expected, we find 
that Hf decreases with decreasing cooling rate. We see that the dependence of Hf on 7 
becomes stronger when we change 7 from small values to intermediate and large values. For 
very large values the dependence becomes, however, weaker again. This can be understood 
by noticing, that for such large values of 7 the cooling process is very similar to a steepest 
descent procedure and that therefore in the limit 7 ^ 00 the dependence of the final state 
on the cooling rate vanishes. 

Also included in Fig. |^ is the result of a cooling process for 7 = 5.0 ■ 10~^ for a system size 
of 2000 particles. This point was generated by averaging over 5 different initial configurations 
and the goal of this "experiment" was to see whether our results are affected in some way 
by finite size effect. Since the result for the system sizes = 1000 and = 2000 are the 
same to within the error bars, we conclude that, for the types of questions investigated in 
this work, finite size effects are not important for the model studied here. 

Since the system will try to decrease its enthalpy during the quench, the cooling process 
can also be viewed as an optimization process in which the system tries to minimize a cost 
function, e.g., the enthalpy. The dependence of the value of the cost function on the amount 
of time spent to find this value has been the subject of quite a few earlier investigations. A 
discussion of some of this previous work can be found in Ref. 0] (see also [| 
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In the just mentioned investigations dealing with coohng or optimization processes it has 
been proposed that the dependence of the cost function, i.e. in our case the enthalpy, on 
the cooling rate is given by a power law or by a logarithmic dependence, i.e. 

Hf{^)=H'f + aa''' (3) 

or 

Hj{^)=H'f + h{-\og^f^ , (4) 

where H^, ai and hi are fit parameters Therefore we have tried to fit our data for 

Hfil) with these functional forms and the best fits are included in Fig. |^ as well. Since 
the functional forms can be expected to hold only for small cooling rates, the three largest 
cooling rates were not included in the fitting procedure. These cooling rates are the ones for 
which the dependence of Hf on 7 becomes weak again, i.e. for which the cooling process is 
similar to a steepest descent. Figure ^ shows that both functional forms are able to fit the 
data very well. Thus within the accuracy of our data and the cooling range investigated it 
is not possible to decide which functional form, if any, is the correct one. 

We also notice that for cooling rates 7 < 2.5 ■ 10"'^ the data points for Hf lie, within the 
accuracy of our data, on a straight line. It is clear that such a functional dependence of Hf 
on 7 cannot be valid for arbitrarily small cooling rates, since this would lead to arbitrarily 
low values of Hf for sufficiently small 7. Thus if there actually exists a regime for which 
there is a linear dependence of Hf on the logarithm of 7, there also must exist at even 
smaller values of 7 a further regime, in which Hf shows the real asymptotic dependence on 
7. However, this regime is currently outside of our computational possibilities. 

In Fig. 1^ we have shown that for large cooling rates the system behaves very anharmonic 
even at low temperatures. Now we find that for these cooling rates the enthalpy of the local 
minimum in which the systems ends at T;, = is relatively large. Thus we conclude that for 
this system the hypersurface of the enthalpy has the property that on average its shape is 
very anharmonic when the enthalpy is large and relatively harmonic when the enthalpy is 
small. 

Similar to the enthalpy we have found that also p/, the density of the system after the 
quench, shows a small but noticeable dependence on the cooling rate. The results are shown 
in Fig. 1^. Note that a change of 7 by about four orders of magnitude leads to a change in 
Pj of only 2%, a figure which is comparable to the one found in experiments (although at 
vastly slower cooling rates) [^. The smallness of the effect shows that it is necessary to vary 
the cooling rate over an extensive range and to use system sizes that are not too small in 
order to obtain results of sufficient statistical accuracy. 

Also in the case of pf we found that, for not too large values of 7, the data can be fitted 
by the two functional forms given by Eqs. @ and (H) equally well. The densities are 1.246 
and 1.253 for the power-law and logarithmic functional forms, respectively. Furthermore we 
see, that, as in the case of the enthalpy, also for this quantity finite size effects are not 
important and that the data points for intermediate and small values of 7 lie on a straight 
line, thus indicating the possibility of a further type of relaxation regime at even lower values 
of 7. 

Since we have now seen that bulk quantities such as the enthalpy and the density show 
a dependence on the cooling rate, it is interesting to investigate how these dependencies 
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can be understood from a microscopic point of view. Therefore we will investigate in the 
following how the microscopic structure of the glass changes when the cooling rate is varied. 
In Fig. 1^ we show the radial distribution functions QAAi"^) sjid gBB{f) for the AA and BB 



correlations |2T[|. The results for gAsif) are similar to the ones for (^aaI^)- The inset in 
Fig. I^a shows that the overall form of gAA{r) does not show a strong dependence on 7. A 
closer investigation of the first nearest neighbor peak showed that for the smallest cooling 
rate its height is about 2% percent higher than the one for the largest cooling rate |T^. A 
somewhat stronger dependence can be observed in the second nearest neighbor peak, part 
of which is shown in the main figure of Fig. ^a. We recognize that the curve for the slowest 
cooling rate (solid bold line) shows more pronounced peaks than the curve for the fastest 
cooling rate (dashed bold line). Thus we conclude that the local order of the system, as 
measured by gAA{f) becomes more pronounced with decreasing cooling rate. 

A similar behavior as the one observed for gAA^r) is found for gBsij) (Fig. |]b). However, 
in this case the height of the first nearest neighbor peak shows a relatively strong dependence 
on the cooling rate. We see that in this case the height of this peak decreases with decreasing 
cooling rate. This behavior can be understood by remembering that esBi the interaction 
energy of the Lennard- Jones potential between two B particles, is only 0.5, as compared 
to 1.5 for eAB- Therefore the system will try to avoid to have two B particles as nearest 
neighbors, and to have A and B particles as nearest neighbors instead, and it will manage 
to do this better the more time it is given to do so, i.e. the smaller the cooling rate is p2[ . 



However, from the figure we also recognize that the height of the first nearest neighbor peak 
decreases when the cooling rate is decreased from very large values to intermediate values, 
but that for even smaller cooling rates the 7-dependence of the height is much weaker, if it 
exists at all. This indicates that the above mentioned driving force for decreasing the peak 
is no longer effective for intermediate and small cooling rates and that therefore a different 
driving force must exist in order to explain the decrease of the enthalpy or the density. 

Having the information on the various radial distribution functions we can now inves- 
tigate the dependence of the different types of coordination numbers of the particles on 
the cooling rate. We define the coordination number 2; of a particle of type a to be the 
number of particles of type (5 that are closer to the first particle than r^f^, the location of 
the first minimum in gaisif)- We have found that the value of r'^^^ depends only we"^'^" — 
the cooling rate |T^ and thus we have chosen the following, 7- independent, values for r. 



r^n = 1-4, r^n = 1-2 and r«« = 1.07. 



a/3 
min' 



In Fig. ^we show the 7 dependence of Pap-, the probability that a particle of type a has 
z nearest neighbors of type (3 for various values of z. We see that the curves for the AA 
pairs. Fig. |^a, for z = 12 and z = 13 increase with decreasing 7, that the curve for z = 10 
shows a decreasing trend and that the curve for 2; = 11 is essentially independent of 7. We 
recognize that the changes that take place in the distribution of the coordination number are 
relatively small, i.e. less than 10%, when the two particles are of type A. This is not the case 
for the BA and BB pairs as is shown by Figs. [pD and c. Here we see that the change in the 
coordination number can be as much as 50% {z = 1 and z = 9 in Fig. ^). Furthermore we 
find that the coordination numbers for BB pairs for z = Q and z = 1 increase and decrease, 
respectively, with decreasing cooling rate (Fig. |^). This observation is in accordance with 
the comments we made in the context of Fig. ^b, where we found that the height of the first 
nearest neighbor peak in the BB correlation function is decreasing with decreasing 7. We 
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also see in Fig. that for values of 7 less than 10"'*, the value of Pbb for 2; = is essentially 
independent of 7, which is in accordance with the comments made before. We also mention 
that from Fig. ^ one should not conclude that, since the curve Paa for z = 12 shows a 
increasing trend, for very small cooling rates the total number of nearest neighbor particles 
of a A particle is 12. What is shown in the figure is just the number of nearest A particles, 
thus the B neighbors are not taken into account. As we will show below the most frequent 
values of the total coordination number for the A particles is 13 and 14. This is different 
from the results one would expect for a one-component system, where at low temperatures 
particles tend to pack locally in an icosahedral structure, thus leading to a coordination 
number of 12. Since in our case we have a mixture of large and small particles it can thus 
be expected that the value of the most frequent coordination number is larger than 12 and 
below we will show that this is indeed the case . 



The fact that the nearest neighbor particles of an A particles do not form an icosahedron 
is also supported by the form of the bond-bond angle distribution function. Here we define 
a bond as the line segment connecting two neighboring particles and a bond-bond angle as 



the angle between two bonds which are connected to a common particle. In Fig. [1^ we show 
this distribution function for the case of three A particles. Note that if the local structure 
of the particles would be a perfect icosahedron we would have only contributions at angles 
63.4° and 116.6°. We see (inset in Fig. p!0|) , that the distribution function shows indeed 
maxima at the mentioned angles, but that the peaks are very broad, thus showing that the 
local structure is not an icosahedron. The overall form of the distribution function shows 
only a small dependence on the cooling rate, the most pronounced effect seems to be that 
the height of the peak near 60° increases and shows the tendency to shift its position to 
larger values of the angle with decreasing cooling rate (main figure). The smallness of these 
effects is in stark contrast with our findings for the network glass former Si02 for which we 
found a strong dependence of the distribution function for the angles on the cooling rate 



Since we have seen now that the arrangement of the particles making up the nearest 
neighbor shell of the A particles is not similar to an icosahedron, we now investigate what 
the nature of this arrangement really is. Furthermore we will also try to understand, what 
the underlying reason for the cooling rate dependence of Hf and pf really is. In order to do 
this we introduce the notion of a "cluster" . We define a cluster of particles as the collection of 
particles given by a central particle and the particles of its nearest neighbor shell (computed 
by using the distances r'^in)- We say that a cluster is of type a^^,^ if the central particle is 
of type a G {A,B} and if it has /i nearest neighbors, u of which are of type B. Note that 
this definition of a cluster can be viewed as a first step of a whole hierarchy of definitions of 
clusters in which, e.g., one defines clusters as the collection of particles that are within the 
first or first two nearest neighbor shell(s), etc. and specifying how many particles of each 
type are in any given shell. However, in order to avoid to complicate the analysis of the 
data too much we restrict ourselves to this simplest kind of cluster definition. 

We define the energy of a cluster to be given by the sum of all pairwise interactions 
between any two members of the cluster and will denote this energy by Ec. The dependence 



of Ec on the cooling rate is shown in Fig. O for various types of clusters. In order not to 



crowd the figure too much, we show only a representative selection of the types of clusters 
we found. (More curves can be found in Fig. 2 of Ref. [^). The dependence of E^. on 7 for 
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the types of clusters not shown is very similar to the one of the types of clusters shown. The 
following observations can be made: 1) For all types of clusters decreases with decreasing 
cooling rate. A variation of 4 decades in 7 gives rise to a decrease of about 1%. 2) For 
a given cooling rate the difference of between the different types of clusters is usually 
much larger than 1%. 3) For a given number of B particles in a cluster, Ec decreases with 
increasing number of A particles in the cluster (as exemplified in the figure by the clusters 
A^^2 (filled symbols)). 4) For a given total number of particles in the cluster, Ec decreases 
with increasing number of B particles (as exemplified in the figure by the clusters A14 (open 
symbols)). Although we have shown here only the dependence of E^ on 7 for clusters of 



type A^ J, similar results hold true also for clusters of the type B^ ,^ p^] . 

Since for each type of cluster a variation of 7 by four decades gives rise to a change 
in Ec on the order of one percent and since the change of the enthalpy, which is related 
to the sum of the energy of the clusters, is also on the order of 1-2%, one might conclude 
that the change of the enthalpy can be rationalized solely by the change of the energy of 
the clusters. Before one can draw this conclusion one has, however, to see whether the 
distribution of the frequency of the different types of clusters is not also changing with the 
cooling rate. In Fig. |12| we show the probability Pa^,^ that a cluster is of type d^^u as a 
function of 7. In order not to crowd the two figures too much, we show only those curves 
for which this probability is not too small. We recognize that there are clusters for which 



this probability is essentially independent of 7, such as, e.g., cluster A132 in Fig. 12a, and 



that there are quite a few clusters for which this probability changes by more than 10% 
(all the clusters that are marked by symbols). Since we have seen from Fig. |lT] that the 
difference between the energy of different clusters is usually much larger than 1-2%, the 
amount of energy a particular cluster changes when 7 is changed by four decades, we thus 
find that the change in the distribution of the frequency of the clusters might be the much 
more important mechanism for lowering the enthalpy than the lowering of the energy of an 
individual cluster. 

Before we address this question further we briefly comment on a different observation. 
From Fig. 12 we recognize, that the observation we made in in Fig. ^a, namely that the curve 
for the AA coordination number 2 = 11 is independent of the cooling rate, does not imply 
that the composition of the clusters containing 11 A atoms is independent of 7. Rather we 
see that the curves for cluster types Ai4^3 and Ai2,i, both of them having 11 A particles, 
show a quite strong dependence on 7. However, this dependence is such that they essentially 
cancel each other, i.e., that the 7 dependence of the abundance of the clusters containing 
11 A particles is almost zero. 

In order to understand the relative importance of the two above mentioned mechanisms 
to lower the enthalpy for the system with decreasing cooling rate, i.e. the lowering of the 
energy of each individual cluster and the changing of the distribution function -P^^^, it is 
useful to look at P{Ec), the probability that a cluster has an energy Ec- In Fig. |13|we show 
this probability for all cooling rates investigated. We see that this distribution shows various 
peaks, each of which can be related to certain types of clusters (see labels). It should be 
noted that the position of a peak is related to the average energy Ec of the corresponding 
cluster and that the height of a peak is related to the frequency with which it occurs. We 
recognize from the figure that the main difference between the distribution for the slowest 
cooling rate (solid bold curve) and the fastest cooling rate (dashed bold curve) is that the 
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height of the various peaks is significantly different. The frequency of cluster types, i.e. 
the value of P{Ec.), with high energies is decreasing and the one of cluster types with low 
energies is increasing with decreasing 7. The location of the peaks, however, shows only 
a very small change in the direction of lower energies, in accordance with point 1) in our 
discussion of Fig. |Tl]. Thus we come to the conclusion that main reason for the observed 
decrease of the enthalpy of the system with decreasing cooling rate is not that the energy of 
the individual clusters is decreasing, but that the distribution of the frequency with which 
these clusters occur shows a relatively strong dependence on the cooling rate. 

Similar to our definition of the energy of a cluster it is also possible to introduce the 
density of a cluster. Therefore it is possible to investigate how the distribution of the 
density of the clusters depend on the cooling rate and we found that, as it is the case for the 
enthalpy, it is the cooling rate dependence of the distribution of the clusters, rather than 
the density of the clusters themselves, that is the main reason for the observed cooling rate 
dependence of the total density of the system |T^. 

To close the subsection containing the results of our analysis of the cooling rate depen- 
dence of the glass at zero temperature, we present the results of our investigation about 
the spectrum of the glass. In order to gain some insight into the nature of the so-called 
boson peak (see, e.g., Ref. ||2^), the low-temperature spectrum of glasses has in recent years 
been the focus of several investigations in which it was attempted to address this question 
by means of computer simulations p6H28|l. Laird and Schober showed that in a glass in 
which the particles interact with a soft sphere potential there exist localized low-frequency 
modes which might be related to the boson peak [^, a result which was also found for a 
model glass for Selenium |2^, and recently Bembenek and Laird have also investigated the 



connection between the shape of the spectrum and the glass transition |28|. In view of this 



work it is therefore interesting to see how the spectrum is affected when the cooling rate 
with which the glass was produced is changed. 

We determined the spectrum of a given relaxed configuration at T = by computing the 
eigenvalues of the dynamical matrix, i.e of d'^V{{ri})/drj^adrk^i3, where j and k are particle 
indices and a and f3 are the cartesian components x, y, z. From each configuration we thus 
obtained 3A^ eigenvalues A, where is the total number of particles, i.e., 1000. From these 
eigenvalues, all of which are positive, since the configuration is locally stable, we computed 
the frequencies u = y^A/m/27r. In Fig. |14| we show the so obtained spectrum Zi^u) for all 
cooling rates investigated. Not included in this figure are the three trivial eigenvalues of zero 
which correspond to a mere translation of the system. A comparison of this figure with the 
spectra shown in Ref. ||2^ for Se shows that the latter spectra are more structured, i.e. show 
more peaks than the one for the Lennard- Jones mixture. This difference is probably due to 
the fact that the potential of the models investigated by the authors of Ref. |^ try to mimic 
directional bonds, whereas the Lennard- Jones model investigated here is more similar to a 
hard sphere model with the corresponding packing structure. This argument is given further 
support by our observation that also the network forming glass Si02 shows a spectrum which 



is much more structured than the one for the Lennard- Jones system presented here ||T3 



We recognize from Fig. |T4| that the spectrum shows only a weak dependence on the 
cooling rate. That there actually is, however, such a dependence is shown in Fig. 15 a, where 
we show an enlargement of the low frequency part of the spectra. We clearly see that with 
decreasing cooling rate the spectrum moves to higher frequencies. A similar behavior is 
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observed for the high frequency wing of the spectrum. Thus we come to the conclusion 
that with decreasing coohng rate the average local environment of the particles transforms 
in such a way that the local potential of the particles becomes stiffer. Such an effect is 
not implausible since we have seen that the density of the system increases with decreasing 
cooling rate and that therefore the particles are closer packed. 

In order to estimate how large of an effect this shifting of the spectrum to higher fre- 
quencies is we proceeded as follows. From Fig. |T3]a we recognize that on the low frequency 
wing of the spectrum the main effect of the cooling rate seems to be that the spectrum is 
shifted to the right, i.e., without a change of the form of the curves. We therefore made 
use of this observation and attempted to shift the curves that have a value of 7 larger than 
3.1 ■ 10~^, the smallest cooling rate, horizontally, i.e., in frequency, by an amount Az/(7) 
such that they collapsed as well as possible with the curve for 7 = 3.1 ■ 10~^ in the range 
where the values of Z^u) is between 0.2 and 1.0. The quality of the collapse between two 
curves was determined with the Spearman rank order correlation coefficient |rT| and the 
shift Ai/(7) was determined by maximizing this coefficient. In Fig. |15|a we show the so ob- 
tained master curve (right set of curves). (In order to make the figure clearer we shifted this 
master curve by 0.02 to higher frequencies.) We thus recognize that it is indeed possible to 
collapse the curves for the different cooling rates onto one master curve, which shows that 
our procedure to estimate the frequency shift of the spectra as a function of the cooling rate 
is reasonable. In Fig. |lpD we show the amount of shifting Az/(7) as a function of 7. We see 
that there is a clear trend that Au decreases with decreasing 7, i.e., that the low frequency 
part of the spectra moves to higher frequencies. It is interesting to note that the curve 
Ai/(7) is approximated quite well by a straight line, i.e., that there is no indication that the 
spectrum converges towards an asymptotic distribution with decreasing 7. This means that 
with respect to this quantity we are not yet seeing the asymptotic behavior expected for 
very small values of 7, thus giving us evidence that there might be, at even smaller values of 
7, a crossover from the 7 dependence of Ai/ we are observing here, to a different one. Note 
that this possibility is also in accordance with our finding in the context of the enthalpy of 
the system at T = 0, where we also had some evidence that there might be a crossover at 
even smaller cooling rates than the ones considered here. 

Although the effect observed here is relatively small, and therefore will probably not 



affect the conclusions drawn in the work of references [B3-E8I, we mention that a similar 



investigation for the strong glass former Si02 has shown that there are systems where the 
spectrum is strongly affected by the cooling rate []T^ . Thus such effects should probably be 
taken into account in the future. 

Apart from this cooling rate dependence of the spectrum we recognize from Fig. |T^ that 
Z^u) also has a system size dependence, since the spectrum shows a gap at low frequencies. 
One knows that acoustic phonons with long wavelength do exist in a glass, as they exist 
in a crystalline solid. Phase space arguments imply that Zlu) oc z/^ for u 0. This is 
not found in the simulation, since the finite linear dimension L of the box implies that the 
largest phonon wavelength is given by L. Hence the corresponding frequency is the lowest 
frequency occurring in the spectrum. Thus this finite size effect should be taken into account 
if one simulates very small systems. 
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IV. SUMMARY AND CONCLUSIONS 



We have performed a molecular dynamics computer simulation of a binary Lennard- Jones 
system in order to investigate the dependence of the properties of a glass on its production 
history. In order to do this we cooled the system, with a cooling rate 7, from its high 
temperature liquid phase to a glass at zero temperature, and investigated how the glass 
transition as well as the so produced glass depend on this cooling rate. 

In qualitative accordance with experiments [0,0 we find that the glass transition becomes 
sharper with decreasing cooling rate, which is also in accordance with simple arguments put 
forward by Angell et al [^]. In particular we see, that in the glass transformation range, 
the specific heat as well as the thermal expansion coefficient show a more rapid change 
the smaller the cooling rate is. Also the glass transition temperature Tg shows a cooling 
rate dependence and we demonstrate that this dependence can be understood by assuming 
that the temperature dependence of the relaxation times of the system is given by either a 
Vogel-Fulcher law or a power-law. 

Our analysis of the glass shows that its density, its enthalpy as well as its structural 
properties, such as the radial distribution function, the coordination numbers and bond-bond 
angles, show a relatively small, but nevertheless noticeable (on the order of a few percent) 
dependence on the cooling rate with which the glass was produced. By investigating the 
cooling rate dependence of nearest neighbor clusters we give evidence that the dependence of 
quantities like the enthalpy or the density is mainly caused by the cooling rate dependence 
of the distribution of the frequency with which the various types of clusters occur. 

Finally we have also investigated how the spectrum Z{v) of the glass depends on the 
cooling rate. We find that Z{v) shows only a weak dependence on 7 and that the main effect 
seems to be that the whole spectrum is shifted towards higher frequencies, a behavior which 
is quite plausible. It is interesting to note that the amount Az/(7) by which the spectrum 
is shifted, does not show any tendency to saturate within the range of cooling rates that we 
are able to access. Therefore it seems that the range of cooling rates falls at least into three 
regimes. The first one includes the very fast cooling rates. In this regime the cooling process 
is essentially a steepest descent procedure and the dependence of the properties of the glass 
on the cooling rate is relatively weak. In the second regime some quantities, such as the glass 
transition temperature or the density, show a relatively strong dependence on the cooling 
rate, but this dependence seems to become weaker with decreasing cooling rate and thus an 
extrapolation to experimental values of the cooling rate seems to be feasible. However, there 
seem to be quantities, such as the shift Az/(7) of the spectrum or, to a lesser extent, the 
enthalpy of the system, which show in this second regime a cooling rate dependence which 
is not possible for arbitrarily small values of 7. Therefore we conclude that there must be a 
third regime of cooling rates in which all quantities actually show a cooling rate dependence 
which is the asymptotic one. 

To conclude we can say that this work has shown that also in computer simulations 
the phenomenon of the glass transition, as well as the properties of the glass, show a clear 
dependence on the cooling rate and that these dependencies can be understood also from a 
microscopic point of view. For the model investigated here these effects are relatively small 
but we know that they might be substantial in other types of systems |]13[. Therefore it 
is advisable that in future simulations which investigate the glass transition, this aspect is 
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not left out since, e.g., the outcome of investigations in which it is tested whether a specific 
model is able to reproduce experimental properties of glasses, might be severely affected by 
the history on how the glass was produced on the computer. 
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FIGURES 



FIG. 1. Enthalpy H of the system versus Tf,, the temperature of the heat bath, for all cooling 
rates investigated, a) Full range of temperature, b) Enlargement of the glass transition region. 
The solid and dashed bold curves are the smallest and largest cooling rates, respectively. 

FIG. 2. Density p of the system versus Tf,, the temperature of the heat bath, for all cooling 
rates investigated. Main figure: The glass transition region. The solid and dashed bold curves are 
the smallest and largest cooling rates, respectively. Inset: Full range of temperature. 

FIG. 3. Specific heat at constant pressure versus Tf, for all cooling rates investigated. The solid 
and dashed bold curves are the smallest and largest cooling rates, respectively, a) Full range of 
temperature, b) Enlargement of the glass transition region. 

FIG. 4. Thermal expansion coefficient Op versus Tf, for all cooling rates investigated. The solid 
and dashed bold curves are the smallest and largest cooling rates, respectively. 

FIG. 5. Glass transition temperature Tg, as determined from the enthalpy (circles) and the 
density (diamonds), versus the cooling rate. The solid and dashed lines are fits with the functional 
forms given in Eq. (|l|) and Eq. (P), respectively. 

FIG. 6. Hf, the value of the enthalpy after the quench, as a function of the cooling rate for 
the system of 1000 particles (open symbols) and the system of 2000 particles (filled symbol). Also 
included are the results of fits with the functional forms given in Eq. (^ (dashed curve) and Eq. (|^ 
(solid curve). 

FIG. 7. pf, the value of the density after the quench, as a function of the cooling rate for 
the system of 1000 particles (open symbols) and the system of 2000 particles (filled symbol). Also 
included are the results of fits with the functional forms given in Eq. d) (dashed curve) and Eq. (0) 
(solid curve). 

FIG. 8. The radial distribution function at T = for all cooling rates investigated. The solid 
and dashed bold line correspond to the slowest and fastest cooling rates, respectively, a) gAA{i")- 
b) 9BB{r)- 

FIG. 9. Probability that a particle has a coordination number z as a function of the cooling 
rate, a) AA pairs, b) BA pairs and c) BB pairs. 

FIG. 10. Distribution of the bond-bond angle between three A particles for all cooling rates 
investigated. Inset: The whole distribution. Main figure: Enlargement of the peak around 60°. 
The solid and dashed bold lines correspond to the fastest and slowest cooling rate, respectively. 

FIG. 11. Ec, the energy of a cluster as a function of the cooling rate for selected cluster types. 
Note that the energy of the cluster of type Afj^^iy decreases with increasing p when v is held fixed 
(filled symbols) and that Ec also decreases with increasing i/ when p is held fixed (open symbols). 
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FIG. 12. Probability to find a cluster of a given type as a function of 7. Only the curves for 
the clusters showing a noticeable change as a function of 7 (marked with symbols) or that are very 
frequent are labeled, a) Clusters around A particles, b) Clusters around B particles. 

FIG. 13. Probability that a cluster has an energy Ec for all cooling rates investigated. The 
solid and dashed bold lines correspond to the smallest and largest cooling rate, respectively. The 
labels of the peaks identify the type of clusters that give rise to the corresponding peak. 

FIG. 14. Z{u), the spectrum of the system at T = for all cooling rates investigated. 

FIG. 15. a) Low frequency part of the spectrum Z{i') for all cooling rates investigated (left set 
of curves). The solid and dashed bold lines correspond to the smallest and largest cooling rate, 
respectively. The right set of curves are the same Z{u) but now shifted by + 0.02 in order 

to collapse them onto a master curve (see text for details), b) The shift as a function of the 

cooling rate. 
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